import numpy
from pyproj import CRS

###################################################
# Geospatial Parameters

# Initiate dictionaries of GEOGRID projections and parameters
#   See http://www.mmm.ucar.edu/wrf/users/docs/user_guide_V3/users_guide_chap3.htm#_Description_of_the_1
projdict = {
    1: 'Lambert Conformal Conic',
    2: 'Polar Stereographic',
    3: 'Mercator',
    6: 'Cylindrical Equidistant'
}
CF_projdict = {
    1: "lambert_conformal_conic",
    2: "polar_stereographic",
    3: "mercator",
    6: "latitude_longitude",
    0: "crs"
}

# Point time-series CF-netCDF file coordinate system
'''Note that the point netCDF files are handled using a separate coordinate system than the grids.
This is because input data are usually in WGS84 or some other global spheroidal datum. We treat
these coordinates as though there is no difference between a sphere and a spheroid with respect
to latitude. Thus, we can choose an output coordinate system for the points, although no
transformation is performed. Properly transforming the points back and forth betwen sphere and
spheroid dramatically increases the runtime of the tools, with clear obvious benefit.'''
pointCF = True  # Switch to turn on CF-netCDF point time-series metadata attributes
pointSR = 4326  # The spatial reference system of the point time-series netCDF files (RouteLink, LAKEPARM). NAD83=4269, WGS84=4326

# Global attributes for altering the sphere radius used in computations. Do not alter sphere_radius for standard WRF-Hydro simulations
sphere_radius = 6370000.0  # Radius of sphere to use (WRF Default = 6370000.0m)
wkt_text = "GEOGCS['GCS_Sphere_CUSTOM',DATUM['D_Sphere',SPHEROID['Sphere',%s,0.0]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]];-400 -400 1000000000;-100000 10000;-100000 10000;8.99462786704589E-09;0.001;0.001;IsHighPrecision" % sphere_radius
#same_spehre = arcpy.SpatialReference(104128)                                   # EMEP sphere (same radius as WRF sphere)

# Unify all coordinate system variables to have the same name ("crs"). Ths makes it easier for WRF-Hydro output routines to identify the variable and transpose it to output files
crsVarname = True  # Switch to make all coordinate system variables = "crs" instead of related to the coordinate system name
crsVar = CF_projdict[
    0]  # Expose this as a global for other functions in other scripts to use
#crsVar = 'ProjectionCoordinateSystem'

###################################################


class WRF_Grid:
    '''
    Class with which to create the WRF-Hydro grid representation. Provide grid
    information to initiate the class, and use getgrid() to generate a grid mesh
    and index information about the intersecting cells.
    Note:  The i,j index begins with (1,1) in the upper-left corner.
    '''
    def __init__(self, rootgrp):
        """
        The first step of the process chain, in which the input NetCDF file gets
        georeferenced and projection files are created
        9/27/2017: This function was altered to use netCDF4-python library instead
        of arcpy.NetCDFFileProperties. Also changed the corner_lat and corner_lon
        index to reflect the lower left corner point rather than lower left center as
        the known point.
        10/6/2017: Proj4 string generation was added, with definitions adapted from
        Adapted from https://github.com/NCAR/wrf-python/blob/develop/src/wrf/projection.py
        5/29/2019: Removed MakeNetCDFRasterLayer_md function call to avoid bugs in ArcGIS
        10.7 (BUG-000122122, BUG-000122125).
        """
        corner_index = 12  # 12 = Lower Left of the Unstaggered grid

        # First step: Import and georeference NetCDF file
        # print('Step 1: NetCDF Conversion initiated...'
        #       )  # Initiate log list for this process

        # Read input WPS GEOGRID file
        # Loop through global variables in NetCDF file to gather projection information
        dimensions = rootgrp.dimensions
        globalAtts = rootgrp.__dict__  # Read all global attributes into a dictionary
        map_pro = globalAtts[
            'MAP_PROJ']  # Find out which projection this GEOGRID file is in
        #print('    Map Projection: %s' % projdict[map_pro])

        # Collect grid corner XY and DX DY for creating ascii raster later
        if 'corner_lats' in globalAtts:
            corner_lat = float(
                globalAtts['corner_lats'][corner_index]
            )  # Note: The values returned are corner points of the mass grid
        if 'corner_lons' in globalAtts:
            corner_lon = float(
                globalAtts['corner_lons'][corner_index]
            )  # Note: The values returned are corner points of the mass grid
        if 'DX' in globalAtts:
            self.DX = float(globalAtts['DX'])
        if 'DY' in globalAtts:
            self.DY = float(globalAtts['DY'])

        # Collect necessary information to put together the projection file
        if 'TRUELAT1' in globalAtts:
            standard_parallel_1 = globalAtts['TRUELAT1']
        if 'TRUELAT2' in globalAtts:
            standard_parallel_2 = globalAtts['TRUELAT2']
        if 'STAND_LON' in globalAtts:
            central_meridian = globalAtts['STAND_LON']
        if 'POLE_LAT' in globalAtts:
            pole_latitude = globalAtts['POLE_LAT']
        if 'POLE_LON' in globalAtts:
            pole_longitude = globalAtts['POLE_LON']
        if 'MOAD_CEN_LAT' in globalAtts:
            #print('    Using MOAD_CEN_LAT for latitude of origin.')
            latitude_of_origin = globalAtts[
                'MOAD_CEN_LAT']  # Added 2/26/2017 by KMS
        elif 'CEN_LAT' in globalAtts:
            #print('    Using CEN_LAT for latitude of origin.')
            latitude_of_origin = globalAtts['CEN_LAT']

        # Check to see if the input netCDF is a WPS-Generated Geogrid file.
        if 'TITLE' in globalAtts and 'GEOGRID' in globalAtts['TITLE']:
            self.isGeogrid = True
        else:
            self.isGeogrid = False
        del globalAtts

        # Handle expected dimension names from either Geogrid or Fulldom
        if 'south_north' in dimensions:
            self.nrows = len(dimensions['south_north'])
        elif 'y' in dimensions:
            self.nrows = len(dimensions['y'])
        if 'west_east' in dimensions:
            self.ncols = len(dimensions['west_east'])
        elif 'x' in dimensions:
            self.ncols = len(dimensions['x'])
        del dimensions

        # Create Projection file with information from NetCDF global attributes
        #sr2 = arcpy.SpatialReference()
        if map_pro == 1:
            # Lambert Conformal Conic
            # if 'standard_parallel_2' in locals():
            #     #projname = 'Lambert_Conformal_Conic_2SP'
            #     #print(
            #         '    Using Standard Parallel 2 in Lambert Conformal Conic map projection.'
            #     )
            # else:
            # According to http://webhelp.esri.com/arcgisdesktop/9.2/index.cfm?TopicName=Lambert_Conformal_Conic
            #projname = 'Lambert_Conformal_Conic_1SP'
            standard_parallel_2 = standard_parallel_1

            projname = 'Lambert_Conformal_Conic'
            Projection_String = ('PROJCS["Sphere_Lambert_Conformal_Conic",'
                                 'GEOGCS["GCS_Sphere",'
                                 'DATUM["D_Sphere",'
                                 'SPHEROID["Sphere",' + str(sphere_radius) +
                                 ',0.0]],'
                                 'PRIMEM["Greenwich",0.0],'
                                 'UNIT["Degree",0.0174532925199433]],'
                                 'PROJECTION["' + projname + '"],'
                                 'PARAMETER["false_easting",0.0],'
                                 'PARAMETER["false_northing",0.0],'
                                 'PARAMETER["central_meridian",' +
                                 str(central_meridian) + '],'
                                 'PARAMETER["standard_parallel_1",' +
                                 str(standard_parallel_1) + '],'
                                 'PARAMETER["standard_parallel_2",' +
                                 str(standard_parallel_2) + '],'
                                 'PARAMETER["latitude_of_origin",' +
                                 str(latitude_of_origin) + '],'
                                 'UNIT["Meter",1.0]]')
            proj4 = (
                "+proj=lcc +units=m +a={} +b={} +lat_1={} +lat_2={} +lat_0={} +lon_0={} +x_0=0 +y_0=0 +k_0=1.0 +nadgrids=@null +wktext  +no_defs "
                .format(str(sphere_radius), str(sphere_radius),
                        str(standard_parallel_1), str(standard_parallel_2),
                        str(latitude_of_origin), str(central_meridian)))

        elif map_pro == 2:
            # Polar Stereographic

            ##            # Set up pole latitude
            ##            phi1 = float(standard_parallel_1)
            ##
            ##            ### Back out the central_scale_factor (minimum scale factor?) using formula below using Snyder 1987 p.157 (USGS Paper 1395)
            ##            ##phi = math.copysign(float(pole_latitude), float(latitude_of_origin))    # Get the sign right for the pole using sign of CEN_LAT (latitude_of_origin)
            ##            ##central_scale_factor = (1 + (math.sin(math.radians(phi1))*math.sin(math.radians(phi))) + (math.cos(math.radians(float(phi1)))*math.cos(math.radians(phi))))/2
            ##
            ##            # Method where central scale factor is k0, Derivation from C. Rollins 2011, equation 1: http://earth-info.nga.mil/GandG/coordsys/polar_stereographic/Polar_Stereo_phi1_from_k0_memo.pdf
            ##            # Using Rollins 2011 to perform central scale factor calculations. For a sphere, the equation collapses to be much  more compact (e=0, k90=1)
            ##            central_scale_factor = (1 + math.sin(math.radians(abs(phi1))))/2        # Equation for k0, assumes k90 = 1, e=0. This is a sphere, so no flattening
            ##
            ##            # Hardcode in the pole as the latitude of origin (correct assumption?) Added 4/4/2017 by KMS
            ##            standard_parallel_1 = pole_latitude
            ##
            ##            printMessages(arcpy, ['      Central Scale Factor: %s' %central_scale_factor])
            ##            Projection_String = ('PROJCS["Sphere_Stereographic",'
            ##                                 'GEOGCS["GCS_Sphere",'
            ##                                 'DATUM["D_Sphere",'
            ##                                 'SPHEROID["Sphere",' + str(sphere_radius) + ',0.0]],'
            ##                                 'PRIMEM["Greenwich",0.0],'
            ##                                 'UNIT["Degree",0.0174532925199433]],'
            ##                                 'PROJECTION["Stereographic"],'
            ##                                 'PARAMETER["False_Easting",0.0],'
            ##                                 'PARAMETER["False_Northing",0.0],'
            ##                                 'PARAMETER["Central_Meridian",' + str(central_meridian) + '],'
            ##                                 'PARAMETER["Scale_Factor",' + str(central_scale_factor) + '],'
            ##                                 'PARAMETER["Latitude_Of_Origin",' + str(standard_parallel_1) + '],'
            ##                                 'UNIT["Meter",1.0]]')

            # 3/30/2020: Testing out utility of providing "Stereographic_North_Pole" or "Stereographic_South_Pole"
            # definitions. However, this may not work for all WRF-supported aspects of stereographic CRS.
            if pole_latitude > 0:
                pole_orientation = 'North'
            else:
                pole_orientation = 'South'
            Projection_String = ('PROJCS["Sphere_Stereographic",'
                                 'GEOGCS["GCS_Sphere",'
                                 'DATUM["D_Sphere",'
                                 'SPHEROID["Sphere",' + str(sphere_radius) +
                                 ',0.0]],'
                                 'PRIMEM["Greenwich",0.0],'
                                 'UNIT["Degree",0.0174532925199433]],'
                                 'PROJECTION["Stereographic_' +
                                 pole_orientation + '_Pole"],'
                                 'PARAMETER["False_Easting",0.0],'
                                 'PARAMETER["False_Northing",0.0],'
                                 'PARAMETER["Central_Meridian",' +
                                 str(central_meridian) + '],'
                                 'PARAMETER["standard_parallel_1",' +
                                 str(standard_parallel_1) + '],'
                                 'UNIT["Meter",1.0]]')

            proj4 = (
                "+proj=stere +units=m +a={} +b={} +lat0={} +lon_0={} +lat_ts={}"
                .format(str(sphere_radius), str(sphere_radius),
                        str(pole_latitude), str(central_meridian),
                        str(standard_parallel_1)))

        elif map_pro == 3:
            # Mercator Projection
            Projection_String = ('PROJCS["Sphere_Mercator",'
                                 'GEOGCS["GCS_Sphere",'
                                 'DATUM["D_Sphere",'
                                 'SPHEROID["Sphere",' + str(sphere_radius) +
                                 ',0.0]],'
                                 'PRIMEM["Greenwich",0.0],'
                                 'UNIT["Degree",0.0174532925199433]],'
                                 'PROJECTION["Mercator"],'
                                 'PARAMETER["False_Easting",0.0],'
                                 'PARAMETER["False_Northing",0.0],'
                                 'PARAMETER["Central_Meridian",' +
                                 str(central_meridian) + '],'
                                 'PARAMETER["Standard_Parallel_1",' +
                                 str(standard_parallel_1) + '],'
                                 'UNIT["Meter",1.0]]')
            proj4 = (
                "+proj=merc +units=m +a={} +b={} +lon_0={} +lat_ts={}".format(
                    str(sphere_radius), str(sphere_radius),
                    str(central_meridian), str(standard_parallel_1)))

        # elif map_pro == 6:
        #     # Cylindrical Equidistant (or Rotated Pole)
        #     if pole_latitude != float(90) or pole_longitude != float(0):
        #         # if pole_latitude, pole_longitude, or stand_lon are changed from thier default values, the pole is 'rotated'.
        #         printMessages(arcpy, ['    Cylindrical Equidistant projection with a rotated pole is not currently supported.'])
        #         ##proj4 = ("+proj=ob_tran +o_proj=latlon +a={} +b={} +to_meter={} +o_lon_p={} +o_lat_p={} +lon_0={}".format(
        #         ##                        sphere_radius,
        #         ##                        sphere_radius,
        #         ##                        math.radians(1),
        #         ##                        180.0 - pole_longitude,
        #         ##                        pole_latitude,
        #         ##                        180.0 + pole_longitude))
        #         sys.exit(1)
        #     else:
        #         # Check units (linear unit not used in this projection).  GCS?
        #         Projection_String = ('PROJCS["Sphere_Equidistant_Cylindrical",'
        #                              'GEOGCS["GCS_Sphere",'
        #                              'DATUM["D_Sphere",'
        #                              'SPHEROID["Sphere",' + str(sphere_radius) + ',0.0]],'
        #                              'PRIMEM["Greenwich",0.0],'
        #                              'UNIT["Degree",0.0174532925199433]],'
        #                              'PROJECTION["Equidistant_Cylindrical"],'
        #                              'PARAMETER["False_Easting",0.0],'
        #                              'PARAMETER["False_Northing",0.0],'
        #                              'PARAMETER["Central_Meridian",' + str(central_meridian) + '],'
        #                              'PARAMETER["Standard_Parallel_1",0.0],'
        #                              'UNIT["Meter",1.0]]')                          # 'UNIT["Degree", 1.0]]') # ?? For lat-lon grid?
        #         proj4 = ("+proj=eqc +units=m +a={} +b={} +lon_0={}".format(str(sphere_radius), str(sphere_radius), str(central_meridian)))

        # # Create a point geometry object from gathered corner point data
        # print('Projection_String:', Projection_String)
        # sr2.loadFromString(Projection_String)
        # sr1 = arcpy.SpatialReference()
        # print('\ncrs_wkt:', wkt_text)
        # sr1 = CRS.from_wkt(wkt_text)
        # sr1.loadFromString(wkt_text)                                                # Load the Sphere datum CRS using WKT
        # point = arcpy.Point()

        # pointX = corner_lon
        # pointY = corner_lat

        # pointGeometry = arcpy.PointGeometry(point, sr1)
        # projpoint = pointGeometry.projectAs(sr2)                                    # Optionally add transformation method:

        # # Get projected X and Y of the point geometry
        # point2 = arcpy.Point(projpoint.firstPoint.X, projpoint.firstPoint.Y)    # Lower Left Corner Point

        # # Populate class attributes
        # self.x00 = point2.X                                                     # X value doesn't change from LLcorner to UL corner
        # self.y00 = point2.Y + float(abs(self.DY)*self.nrows)                    # Adjust Y from LL to UL
        # self.proj = sr2
        # #self.WKT = Projection_String.replace("'", '"')                          # Replace ' with " so Esri can read the PE String properly when running NetCDFtoRaster
        # self.WKT = sr2.exportToString()
        # self.proj4 = proj4
        # self.point = point2
        # self.map_pro = map_pro
        # printMessages(arcpy, ['    Georeferencing step completed without error.'])
        # del point, sr1, projpoint, pointGeometry

        self.Projection_String = Projection_String

    def GeoTransform(self):
        '''
        Return the affine transformation for this grid. Assumes a 0 rotation grid.
        (top left x, w-e resolution, 0=North up, top left y, 0 = North up, n-s pixel resolution (negative value))
        '''
        return (self.x00, self.DX, 0, self.y00, 0, -self.DY)

    def GeoTransformStr(self):
        return ' '.join([str(item) for item in self.GeoTransform()])

    def getxy(self):
        """
        This function will use the affine transformation (GeoTransform) to produce an
        array of X and Y 1D arrays. Note that the GDAL affine transformation provides
        the grid cell coordinates from the upper left corner. This is typical in GIS
        applications. However, WRF uses a south_north ordering, where the arrays are
        written from the bottom to the top.
        The input raster object will be used as a template for the output rasters.
        """
        print('    Starting Process: Building to XMap/YMap')

        # Build i,j arrays
        j = numpy.arange(self.nrows) + float(
            0.5)  # Add 0.5 to estimate coordinate of grid cell centers
        i = numpy.arange(self.ncols) + float(
            0.5)  # Add 0.5 to estimate coordinate of grid cell centers

        # col, row to x, y   From https://www.perrygeo.com/python-affine-transforms.html
        x = (i * self.DX) + self.x00
        y = (j * -self.DY) + self.y00
        del i, j

        # Create 2D arrays from 1D
        xmap = numpy.repeat(x[numpy.newaxis, :], y.shape, 0)
        ymap = numpy.repeat(y[:, numpy.newaxis], x.shape, 1)
        del x, y
        print(
            '    Conversion of input raster to XMap/YMap completed without error.'
        )
        return xmap, ymap

    def grid_extent(self):
        '''
        Return the grid bounding extent [xMin, yMin, xMax, yMax]
        '''
        xMax = self.x00 + (float(self.ncols) * self.DX)
        yMin = self.y00 + (float(self.nrows) * -self.DY)
        #extent_obj = arcpy.Extent(self.x00, yMin, xMax, self.y00)
        return (self.x00, yMin, xMax, self.y00)

    def xy_to_grid_ij(self, x, y):
        '''
        This function converts a coordinate in (x,y) to the correct row and column
        on a grid. Code from: https://www.perrygeo.com/python-affine-transforms.html
        Grid indices are 0-based.
        '''
        # x,y to col,row.
        col = int((x - self.x00) / self.DX)
        row = abs(int((y - self.y00) / self.DY))
        return row, col

    def grid_ij_to_xy(self, col, row):
        '''
        This function converts a 2D grid index (i,j) the grid cell center coordinate
        (x,y) in the grid coordinate system.
        Code from: https://www.perrygeo.com/python-affine-transforms.html
        Grid indices are 0-based.
        '''
        # col, row to x, y
        x = (col * self.DX) + self.x00 + self.DX / 2.0
        y = (row * -self.DY) + self.y00 - self.DY / 2.0
        return x, y


def add_CRS_var(rootgrp,
                sr,
                map_pro,
                CoordSysVarName,
                grid_mapping,
                PE_string,
                GeoTransformStr=None):
    '''
    10/13/2017 (KMS):
        This function was added to generalize the creating of a CF-compliant
        coordinate reference system variable. This was modularized in order to
        create CRS variables for both gridded and point time-series CF-netCDF
        files.
    '''

    # Scalar projection variable - http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/reference/StandardCoordinateTransforms.html
    proj_var = rootgrp.createVariable(CoordSysVarName,
                                      'S1')  # (Scalar Char variable)
    proj_var.transform_name = grid_mapping  # grid_mapping. grid_mapping_name is an alias for this
    proj_var.grid_mapping_name = grid_mapping  # for CF compatibility
    proj_var.esri_pe_string = PE_string  # For ArcGIS. Not required if esri_pe_string exists in the 2D variable attributes
    proj_var.spatial_ref = PE_string  # For GDAl
    proj_var.long_name = "CRS definition"  # Added 10/13/2017 by KMS to match GDAL format
    proj_var.longitude_of_prime_meridian = 0.0  # Added 10/13/2017 by KMS to match GDAL format
    if GeoTransformStr is not None:
        proj_var.GeoTransform = GeoTransformStr  # For GDAl - GeoTransform array

    # Projection specific parameters - http://www.unidata.ucar.edu/software/thredds/current/netcdf-java/reference/StandardCoordinateTransforms.html
    if map_pro == 1:
        # Lambert Conformal Conic

        # Required transform variables
        proj_var._CoordinateAxes = 'y x'  # Coordinate systems variables always have a _CoordinateAxes attribute, optional for dealing with implicit coordinate systems
        proj_var._CoordinateTransformType = "Projection"
        proj_var.standard_parallel = sr.standardParallel1, sr.standardParallel2  # Double
        proj_var.longitude_of_central_meridian = float(
            sr.centralMeridianInDegrees
        )  # Double. Necessary in combination with longitude_of_prime_meridian?
        proj_var.latitude_of_projection_origin = float(
            sr.latitudeOfOrigin)  # Double

        # Optional tansform variable attributes
        proj_var.false_easting = float(
            sr.falseEasting
        )  # Double  Always in the units of the x and y projection coordinates
        proj_var.false_northing = float(
            sr.falseNorthing
        )  # Double  Always in the units of the x and y projection coordinates
        proj_var.earth_radius = sphere_radius  # OPTIONAL. Parameter not read by Esri. Default CF sphere: 6371.229 km.
        proj_var.semi_major_axis = sphere_radius  # Added 10/13/2017 by KMS to match GDAL format
        proj_var.inverse_flattening = float(
            0
        )  # Added 10/13/2017 by KMS to match GDAL format: Double - optional Lambert Conformal Conic parameter

    elif map_pro == 2:
        # Polar Stereographic

        # Required transform variables
        proj_var._CoordinateAxes = 'y x'  # Coordinate systems variables always have a _CoordinateAxes attribute, optional for dealing with implicit coordinate systems
        proj_var._CoordinateTransformType = "Projection"
        proj_var.longitude_of_projection_origin = float(
            sr.longitudeOfOrigin
        )  # Double - proj_var.straight_vertical_longitude_from_pole = ''
        proj_var.latitude_of_projection_origin = float(
            sr.latitudeOfOrigin)  # Double
        proj_var.scale_factor_at_projection_origin = float(
            sr.scaleFactor)  # Double

        # Optional tansform variable attributes
        proj_var.false_easting = float(
            sr.falseEasting
        )  # Double  Always in the units of the x and y projection coordinates
        proj_var.false_northing = float(
            sr.falseNorthing
        )  # Double  Always in the units of the x and y projection coordinates
        proj_var.earth_radius = sphere_radius  # OPTIONAL. Parameter not read by Esri. Default CF sphere: 6371.229 km.
        proj_var.semi_major_axis = sphere_radius  # Added 10/13/2017 by KMS to match GDAL format
        proj_var.inverse_flattening = float(
            0
        )  # Added 10/13/2017 by KMS to match GDAL format: Double - optional Lambert Conformal Conic parameter

    elif map_pro == 3:
        # Mercator

        # Required transform variables
        proj_var._CoordinateAxes = 'y x'  # Coordinate systems variables always have a _CoordinateAxes attribute, optional for dealing with implicit coordinate systems
        proj_var._CoordinateTransformType = "Projection"
        proj_var.longitude_of_projection_origin = float(
            sr.longitudeOfOrigin)  # Double
        proj_var.latitude_of_projection_origin = float(
            sr.latitudeOfOrigin)  # Double
        proj_var.standard_parallel = float(sr.standardParallel1)  # Double
        proj_var.earth_radius = sphere_radius  # OPTIONAL. Parameter not read by Esri. Default CF sphere: 6371.229 km.
        proj_var.semi_major_axis = sphere_radius  # Added 10/13/2017 by KMS to match GDAL format
        proj_var.inverse_flattening = float(
            0
        )  # Added 10/13/2017 by KMS to match GDAL format: Double - optional Lambert Conformal Conic parameter

    elif map_pro == 6:
        # Cylindrical Equidistant or rotated pole

        #http://cfconventions.org/Data/cf-conventions/cf-conventions-1.6/build/cf-conventions.html#appendix-grid-mappings
        # Required transform variables
        proj_var.grid_mapping_name = "latitude_longitude"  # or "rotated_latitude_longitude"

        # No extra parameters needed for latitude_longitude

        # If rotated pole:
        #printMessages(arcpy, ['        Cylindrical Equidistant projection not supported.'])
        #raise SystemExit

    # Added 10/13/2017 by KMS to accomodate alternate datums
    elif map_pro == 0:
        proj_var._CoordinateAxes = 'lat lon'
        proj_var.semi_major_axis = sr.semiMajorAxis
        proj_var.semi_minor_axis = sr.semiMinorAxis
        if sr.flattening != 0:
            proj_var.inverse_flattening = float(
                float(1) / sr.flattening)  # This avoids a division by 0 error
        else:
            proj_var.inverse_flattening = float(0)

    # Global attributes related to CF-netCDF
    rootgrp.Conventions = CFConv  # Maybe 1.0 is enough?
    return rootgrp
